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^1^, Abstract 



< 



It is often not possible to construct a probability density function that describes the data. This can happen if there 
^SJ ' is no analytic description, and the number of parameters is too large so that it is impossible to simulate and tabulate all 

combinations. In these situations it is still interesting to rank simulation sets performed with different parameters in how 
well they compare to data. We propose a solution that appears to be better suited to this task than some of the obvious 
alternatives. 
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It is often the case that the mean rate of counts in a bin is not known exactly but rather approximated with simulation. The 
simulation can be repeated many times, obtaining a total number of counts of s in Ug trials, and the expected rate of counts 

C3 ' is often approximated as /i = s/ng. For the sake of generality let's assume that we repeat the experiment n^ times and 
collect a total of d counts. In order to fit for some unknown property of the experiment one often maximizes the likelihood, 
and for convenience that is usually done by minimizing the minus log likelihood, — In C The minus log likelihood based 

^ , on the Poisson probability of our observation is given by 

in 

CO ' — In £ = In d! + UdiJ- — d ■ hi{ndfJ.) ■ 

i>: 

^^ ■ Unfortunately this expression is only an approximation as the quantity /i is not known precisely but was calculated from 
■^ I simulation and is known within statistical uncertainties corresponding to the total number of simulated counts s in our bin. 
^D ■ When the counts s in simulation and d in data are large, one can write an expression similar to the x^'- 
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s/n1 + d/n^ 



where the total uncertainty in the denominator is computed as the square root of the sum of squares of mutually independent 
j^ \ statistical uncertainties of s/ug and d/ud- 

One may approach this problem from the Bayesian point of view: the counts in simulation are distributed with a Poisson 
probability around some unknown value of the true rate /i, so we convolve that probability (treating it as the likelihood) 
with the probability to observe the d counts in data, and with some prior (taken as /z^ in the following expression): 

-ln£ = -ln^ ;.d^ 

-In^ ; —- -din = sin h (z + l)ln{ns +nd). 

s\ ■ d\ Us + Ud ris + Ud 

A choice that is often used is z = — 1, corresponding to a non-normalizable prior of /i^^. 

In the next section we introduce a possible new treatment that considers statistical uncertainties in both data and simu- 
lation, and appears to perform better than the alternatives listed above, as demonstrated later in section |9] 



2 Likelihood description with statistical uncertainties only 

Consider a repeatable experiment that is performed n^ times to collect a total of d counts with a per-event expectation of Hd 
(we call a single instance of this experiment an "event"). We predict the result of the experiment with the simulation, which 
collects s counts in Us simulated events and a per-event expectation of /i,,. 

Given that the total count in the combined set of simulation and data is s + d, the conditional probability distribution 
function of observing s simulation and d data counts is 






An obvious constraint that is implied here is ns^g+n^^d ~ s + d, which can be derived from the normalization requirement 

s+d 
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If the data and simulation are completely unrelated the best possible estimates of ^s and ^d are determined by maxi- 
mization of the probability function given above with the constraint Usiis + ridUd ^ s + d, which results in the estimates 

s d 

Ms = — : Md = — ■ 
ris Ud 

Now, the alternative hypothesis that we could try to verify is that data and simulation counts are described by the same 
process, i.e., come with the same mean ji = ^g = Md- This identity together with the constraint ngfis + «dMd — s + d 
uniquely determines the values of 

s + d 

^j. = ^is = fJ-d = ■ • 

Us +nd 

We can now compare the probabilities of the two of the above hypotheses by forming a likelihood ratio 

P(same process) f ris / ^ \ f iT-d / d \ f 1^ \ f ^^ \ 

P (independent processes) \ns + Ud s + d) yug + Ud s + d J \s/nsj \d/ndj 

This ratio is the "likelihood" that should be used for comparison of different simulation sets with data. The denominator 
in the above expression effectively factors out the dependence of the probability P on the number of "states" (i.e., possible 
values of s) around the most likely value of s. (This dependence comes in because of the normalization condition: the sum 
of P over all states should be 1, which means that the value of P for the most likely values of s is lower for larger s + d.) 

3 Generalization to many bins 

If there are several bins {i} in which simulation and data counts are compared, the conditional probability can be written as 

{s + D)\ rr ^ "«^n '' TT ^ "rf/^d ""* 



.(Ki, Ml: (.1. ims . m ^ nj^fi^ . n ^i^j . n ^1^ 

Here we use notations S = X^i ^i, D ~ '^- di. The probability sum of 1 requires '^iirisiJ-l + nd^^d) = S + D. Taking 
the negative logarithm, losing constant terms, and introducing a Lagrange multiplier term for this constraint (with a new 
unknown Q, this becomes: 

F = - ^ s, ln(M>.) -^d, Hiilind) + C ' (^ "''^^^ + J2 "^dfJ-d - S - D). 

i i i i 

If data and simulation are independent, this expression is minimized for each ^\, /i^ independently: 

dF Si ^ dF di ^ ■ Si : di 

'+C-"s = 0, TTT^ - + C,-nd = Q +> M. = 7 ' ^'d 



di^l Ml' ' dfi^, M^i ' "^ C-n,' "" C-fid- 



Plugging this back into the constraint equation (which we also get back by setting dF/dC, = 0), we get 

If data and simulation come from the same distribution we require /i' = /i* = /i^, and minimize against /x*: 

dF Si + di , . , , . „ i s, + d, 

+ C ■ ("s + rid) = ^ ^' 



Once again, plugging this back into the constraint relation we get 

Thus, the expressions derived in the previous section for 1-bin situation are valid per-bin when there are more than one 
bin, and we get back the likelihood ratio formula 

P(same process) TJ f ^^^ \ ' TJ f ^^ \ ' ■ u « Si + di 

P(independent processes) ^^\si/nsj \di/nd) ' ng + Ud 

We compare the performance of reconstruction using this formula with that using UkeUhood expressions listed in the intro- 
duction in section |9] 

4 Likelihood description: adding model errors 

The error in describing data with simulation (i.e., describing fjid with /is) is often non-zero. In such a case one may quantify 
the amount of disagreement between data and simulation with a x^: 

X — _i 



Instead of setting /i^ = /i^ as in the previous sections we assume that a difference between /i^ and /i^ can exist due to this 
systematic error and is modeled with a likelihood penalty term 

In^ (/id/Ms) 

The likelihood ratio is therefore determined as 

P(same process) _ ( ^s \'' f fJ-d \ ln^(/^d/Ms) 



P(independent processes) \s/nsj \d/ndj ~2a- 

where the /i^ and ^d are determined by maximizing 

with the constraint UgHs + ndUd — s + d. Taking the negative logarithm, losing constant terms, and introducing a Lagrange 
multiplier term for this constraint (with a new unknown Q, this becomes: 

-sln(/isns) - d\n{fj,dnd) + 7-^ In h C ' ("sMs + nd^id - s - d) = F. 

2a^ fig 

The function F{fis, IJ-d) can be easily minimized against /i^ and /i^, yielding estimates of these quantities. To demon- 
strate this, first the derivatives of F are calculated and set to 0: 

OF 1 , Aid n 



MrfTj — = Qt^-dTid - a! + ^j In — = 0. 

The sum of these, Q ■ {fisTis + jJ^dnd) ~ s + d, results in the value for C = 1- The derivative of F with respect to C gives 
back the constraint Usfis + nd^J■d ~ s + d, which yields an expression of fid as a function of fig- Plugging it into the first of 
the above two equations one gets 

/ = Mst; — (A^s, MdlA^sj) = Ms"-s - s ^ In = 0. 

Ofls <J Us 

This equation can be solved with a few iterations of the Newton's root finding method starting with a solution to 

s + d 



fJ-s ^ fJ-difJ-s)- fJ-s = fJ-d 



ris + Ud 
At each iteration the value of Hs is adjusted by —f/f, where the derivative is evaluated as 






Once the likelihood function is solved for the best values of ^s and fid, these can be plugged into the likelihood ratio 
given above. One can now write the likelihood ratio as a sum over all bins: 



ln^ = E 



,anfl/^ + d,ln^i/^ + ^ln^^ 
fll /4 2cr^ fls 



This is an improved expression compared to the one used in |l3]. The probability P(same process) can also be thought of 
as a convolution of the binomial probability part of the expression with the penalty term. Solving the convolution integral 
approximately with the Laplace's method results (up to a constant term) in an expression for P(same process) given above. 

5 Likelihood description of data with weighted simulation 

One can apply the method for calculating the likelihood ratio of the previous section to a situation that is common when the 
number of data counts dk in bin k measured during time td is fitted with a number of simulation counts Sfc^ , each representing 
a possibly different (for weighted simulation) time tki (usually related to the event weight Wki as Wki ■ tki — td). Although 
we can assume that all Ski = 1 without the loss of generality, we continue with the notation Ski- The combined number 
of events in data and simulation is then S + D, where S = J2k ^k, Sk = J2i ^ki, D = ^j, dk- The expression for the 
conditional probability is now 

P{{fif}A^4}\{sk^}Adk}\s + D 

The probability sum of 1 necessitates the constraint 



ki k 

Taking the negative logarithm of P, losing constant terms, and introducing a Lagrange multiplier term for this constraint 
(with a new unknown Q, we get: 



^ T\k^^k^-T{kdk\ \\\s + d) 


TT ( tdt^d 
\^{S + D 


constraint 




Y,tk^^J^'T + Y.^d^i''d = s + D. 





ki k \ ki k 
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If data and simulation are independent, this expression can be minimized for each /i**, fi^ independently: 

OF ''k^ dF _ dk k^ _ Sk^ fc _ dk 



Plugging this back into the constraint equation (which we also get back by setting dF/dC, = 0), we get 



E'«7^ + E'^7^-Ei^-« + ^ - «-!■ 



If data and simulation come from the same distribution we require /i^ — ^^ /x^' for each k. These conditions can be 
introduced into the above expression for F as additional terms (with new unknowns ^fc): 

ki k \ ki k / k \ i / 

Derivatives with respect to ( and ^fe give back the constraint equations. The other derivatives are: 

dF Ski ... „ dF dk 
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Multiplying the first equation by /i^*, the second by /i^, and summing them together, we get 

ki k \ ki k / k \ i / k 

Therefore (^ = 1. To find ^ we substitute the expressions for /i^* and /ijj into constraints for ^fc: 

Ski k ^k , V^ Ski dk 



M 



ki 



Md 



tki + £,k 

Therefore, the likelihood ratio is 

P(same process) 



td — £,k 



E 



tki + S,k td — ^k 



P(indcpendent processes) 
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tki 



tki + (,k 



n 



td 



td — £,k 



d,. 



The equation for ^^ is similar to equation 15 of H]. As suggested there, we solve them for each k by Newton's method 
starting with ^^ = 0, ensuring that - mm[s^^>o}itki) < £,k < td'. 



fk = 1/ 
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Ski 



tki + S,k 



td — Cfc dfk _ MT-^ 

dk ' d^'^ 



Ski 



{tk^+ik)' 



r/ 
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Ski 



tki + ^k 
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dk 



^fe(next) =^fe- 



fk 



dfk/dik 



The particular form of function fk above (inverted compared to the original equation for ^k) was chosen to linearize the 

problem in simple cases (e.g., all simulated events having the same weight). After the first iteration (starting with f^ = 0) 

we get 

. td/dk - 1/I]i Ski/tki 1 - dk/ruk 

— td 



l/dk + J:^Sk^/tl/[J:.,Sk^/tk^? ^" l + dk-sl/ml 

with nik ^^^SkiWki, ^l^^^Skiwl^, Wik 



ik 



Here mk is the total simulation in bin k evaluated for time td (i.e., is a prediction of data in bin k)\ Sk is a statistical 
uncertainty of the value nik- If irik and dk are not too far from each other, and the statistical uncertainty of the simulation is 
much smaller than that of data, i.e., Sk <C \/dk, then Cfe ~ id • (1 — dk/nik)- This holds in the limit of infinite simulation 
statistics, when all tki —> oo, and thus, Wki -^ 0. Continuing with this approximation. 



P(same process) 
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P (independent processes) -'■-'- \ 1 + Wki • (1 — dk/mk) 



n T^ 
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(1 - dk/mk) 



exp 



- ^ SkiWki ■ (1 - dk/mk) 

ki 



exp(dfc - mk) 



mk 
dk 



dk 



Up to a constant term this is a product of the usual expressions for the Poisson UkeUhood, 

p{m; d) = —— ■ exp(— m). 

However, the expression above is also precisely identical to the product of ratios p(to; d)/p{d; d), where m = d in the 
denominator is the value at which the Poisson likelihood p(to; d) achieves its maximum. 

We would like to emphasize that using the usual Poisson expression for the likelihood fit of data to simulation is, strictly 
speaking, only correct in the limit of counts of the simulation being infinite in all of the bins used in the likelihood fit, and 
all of the simulation event weights being infinitesimally small. An uneven distribution of simulation counts in bins (e.g., 
energy bins) may create a bias in the resulting fit (i.e., bias in energy). In this section we advocate using the exact expression 
for the likelihood ratio given above. We are aware of several recent works (e.g., [4|) trying the new expression and finding 
that it works better than the usual Poisson Ukelihood. 

6 Likelihood description with weighted simulation and model errors 

Instead of assuming the condition of '^- /i^* = /i^ used in section|5]we add a penalty term similar to section|4l 

To simplify the following expressions we introduce notations /i^ = J^i I^T ^i^^ ^k = ln(Ms//^d)- '^°' i^ '■^^ '^^'■^ ^'^'^ 
simulation come from the similar distributions, discrepancies being described by the above penalty term, the expression for 
Fis: 

ki k \ ki k / k 

The derivatives with respect to the unknown /i^', /ijj are: 

dF Ski , /- ^ , Cfe n ^^ d,k , ^ _, Ck 



+ C-t.. + -^ = o, _ = -^ + c.t,--^ = o. 



As in section ID we prove that C = 1 by multiplying these equations by /x^' and /i§ and summing them together. We can 
rewrite the above equations to get expressions for /i^' and /i^: 

ki Skz k ^ f^ ,^k 

tk^ + ^ td\ 

If we find £^k we could determine all of /i^'* and fi'^: first, /i§ is given by the second of the above equations, then we use 
A*s = ^fcMd ^nd determine /i^* from the first of the above equations. We can now write an equation for ^k'- 

k V^ ki V^ ^ki V^ Ski -, , k ^ I J , ik 

Similar to section |5] this equation can be solved with Newton's method starting with ^^ = 1. We can now write the 
hkelihood ratio as 

P (similar process) yj fhtfJ-s^V"^ Tf f^df^d\ ' TT f '^ c2 



P (independent processes) \^ \ Ski J ^}- \ dk J \^ \ Scr-^ 
To find optimal parameters of the simulation one typically maximizes this by minimizing the minus log likelihood: 

We note that the derivatives of this expression with respect to parameters tki are equal to the partial derivatives of this 
expression with respect to parameters tki, since the rest has terms proportional to the derivatives with respect to /i^* and /ijj, 
which were set to by the above effort. This could be useful when searching for the optimal values of tki themselves. 
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Figure 1 : Dependence of the one-term likelihood on the simulation weight. The data and simulation are both set to 1 count, 
so the minimum is at w = 1, as expected. Without noise the likelihood expressions of both sections|5](o' = 0) and|6]((T = 1) 
diverge at (i.e., at large negative values of logj^Q w). Adding noise (with rate 77 = 10~^) results in a finite value at for 
both likelihood expressions. Also note that the curves with a = 1 have wider minima, as expected. 



7 Adding noise 

The likelihood expressions of the previous two sections can only be computed when all bins that have data counts contain 
at least one or more simulated counts. If there are bins with data but no simulation there must have been some effect that 
was not simulated, or the amount of simulation was insufficient. In such cases the likelihood calculation is impossible 
unless one excludes these bins from the consideration. For a meaningful comparison of likelihood values computed for 
different simulation sets all bins that have no simulated counts in at least one of the simulated sets must be excluded from 
all likelihood evaluations. 

It is sometimes desired to not simulate a simple effect of a known rate rjk (i.e., noise) but rather directly introduce j]k into 
the likelihood calculation. This would coincidentally solve the problem with bins that have no simulated counts described 
in the preceding paragraph. The constraint for rates of section|5]is modified as J^i Ms* + Vk — M^' ^'^'^ '^^e penalty term of 
sectionals modified via a new expression for ^fe: ^^ — ln[(/i^ + '7fc)/Md]- When the value rjk of noise rate is known much 
better than the a or the penalty term assumes, the new expression for ^fc may be written as ^fe = ln[/i^/(/i^ — r/k)] instead. 
Although the expression for (^ becomes more complicated, it can be still approximated as 1 in the case of many bins with a 
comparably large total number of counts S + D. Figure[T]compares the behavior of the likelihood functions (with only one 
simulation and one data counts) of sections|5]and|6]with and without noise. 



8 Binning 

If the data is binned in equal-length consecutive bins one may occasionally benefit from combining the counts in contiguous 
subsets of bins into larger bins of length L (L is the number of initial-size small bins in the larger bin). The benefit of 
combining the counts of several bins must be weighted against the possibility of measuring the change in rate between the 
bins. If the counts are comparable to each other or are very small in the two adjacent bins (i.e., are within the statistical 
uncertainties), then the rate change across the these bins is consistent with zero, and thus the bins should be combined. This 
is the principal idea behind the "Bayesian blocks" algorithm of ^. Below we summarize the algorithm, as it is used in 
section |9] Note how the probability ratio construction used by this algorithm is quite similar to the probability ratio used in 
this paper 



The multinomial probability for a particular distribution of counts in bins under consideration is 



D\ 






y^^^ 



The first hypothesis is that the counts in these bin are independent from each other It follows that jii — di for all i. The 
second hypothesis is that the counts in all of the bins under consideration are governed by the process with the same rate 
II = Hi for all i. If follows that /u = D/L, where L is the number of bins. 

To determine whether the bins should be combined or be considered separately we form the likelihood ratio: 

P(onebin) (D/L)'^ , v-^ , 



di ' ""^^^ ^ /_^' 



P(many bins) Yl d 

Now we note that this ratio never exceeds 1, as the value of the likelihood maximized with a constraint /i = fM is 
necessarily smaller or equal to than the value of the likelihood maximized without any restrictions. 

Note that for any permutation of bins i the ratio above does not change. However, with the assumption of the same 
rate across the bins all (LI) of the permutations are equivalent to merely being different statistical realizations of the same 
experiment, all leading to the same term Y[i d.i ' in the denominator. On the other hand there is only one possible realization 
of the data arrangement that leads to the term (D/L)^ in the nominator 

We will now make an intuitive conclusion that it is fair to decide whether to combine the bins i into a single large bin 
with length L based on a comparison of the ratio above with 1/Ll. We can easily generalize this rule for comparing various 
arrangements of bins k of different sizes Lk'- we choose the arrangement that leads to the largest sum (taking the log of the 
probability ratio) 

This is different from the consideration given in [2] in the expression for the last term, which in 1*21 is given by « — ln(8). 
This reference also contains a description of the algorithm for the search of the best bin arrangement, which is easy to 
implement and is very efficient. 

To get a little more insight into the expression given above we write it for a number L of bins i each with approximately 
the same counts di — d + Si, D — J^i di = d ■ L, 6i/d <C 1: 



E^ 



^d, •lnd, + ln(l!) = ^(d + J,)ln(d + ^.)«^fdlnd+|,) = ^^'^ (t) + IT.^ 



We compare this with the expression for a single large bin of length L: 

D\n(^] +ln(L!). 



So, the decision to combine the bins is based on whether 



IE 



X2 / T \ s:2 

" < ln(L!) « L • In - - L, or ^ -^ < ^, 



d 

i.e. when the x^ (i.e., goodness-of-fit) of description of data with the mean d in the bins is on the order of or smaller than 
the number of bins (i.e., degrees of freedom of the x^). 

9 Comparison with other methods 

In this section we will compare the performance of parameter reconstruction using the likelihood function introduced in 
section[3]and three likelihood approaches enumerated in the introduction: Poisson, x^, and Bayesian. We also reconstruct 
using the "true" Poisson likelihood by calculating the expectation in each bin directly from the probability distribution 
function (PDF). We will take the distribution of counts following the following PDF: exp(— a;//i)//i with fio = 5. The sets 
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Figure 2: Left: data and simulation for Ud = 100 and Us = 100. Also shown in blue solid line is the binning optimized 
with the method described in section [8] Right: likelihood values for random simulations in the vicinity of the true value of 
Mo = 5. 



are drawn with the mean total number of counts per set of 10. This "experiment" is performed n^ = 1, 10, 100, or 1000 
times, and each simulation is performed Ug — 1, 10, 100, 1000 times. 

The data and simulation are optionally rebinned according to the method described in section [8] We investigated the 
binning strategy optimized for combined counts of data and simulation, which often gave the best result, and compared it 
to the binning optimized for counts of just data, which is often almost as good but more practical (as it is only done once 
per data sample). The reconstruction of the value of parameter fi is performed by simulating at points randomly and log- 
uniformly sampled in the range /io/5 to 2/io and comparing resulting likelihood values. A typical result of such sampling 
procedure is shown in Figured along with a comparison between a representative data sample and a simulation sample at 
the best fit point. 
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Figure 3: Distributions of the reconstructed value of /i for 1000 random drawings of data with n^ = 100. For each data set, 
200 simulation sets with /i sampled log-uniformly in the interval (/io/5, 2/io) ^e drawn with Ug — 10 and the one with the 
best likelihood value is chosen to represent the result of the fit. The top row uses the bins of size 1, the bottom row uses 
the bins with size optimized according to the method of section|8] The Poisson "true PDF" likelihood is always the best, as 
expected, and does not depend on the binning method. In all other cases the result is improved with the optimized binning. 
The "full likelihood" corresponds to the likelihood description of this paper and gives the best result in this example. 



We compare our findings in Table [T] A representative set of distributions of the reconstructed values of /i is shown in 
Figure |3]for s = 10, and d = 100. The new likelihood gives the best result (closest to fiQ = 5 with smallest spread) more 
often than the other methods, except, of course, the method using the Poisson likelihood with the mean estimated from the 
analytical expression for the PDF. We also comment that the optimized binning improves the result in most cases. 
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4.86 ± 0.65 


3.35 ± 0.51 


4.41 ± 0.92 


4.05 ± 0.89 


5.01 ±0.52 


100 


100 


4.94 ± 0.35 


4.95 ± 0.34 


4.63 ± 0.34 


4.72 ± 0.36 


5.00 ±0.16 


100 


1000 


4.93 ± 0.26 


5.39 ± 0.33 


4.71 ± 0.26 


4.86 ± 0.28 


5.00 ± 0.06 


1000 


1 


4.87 ± 1.66 


1.44 ± 1.09 


6.88 ±3.21 


3.11 ±1.77 


4.94 ±1.67 


1000 


10 


4.94 ± 0.56 


2.25 ± 0.37 


4.51 ±0.82 


4.72 ± 0.70 


5.01 ±0.50 


1000 


100 


4.98 ± 0.23 


4.46 ± 0.23 


4.76 ±0.23 


4.93 ± 0.24 


4.99 ±0.17 


1000 


1000 


4.99 ± 0.13 


4.99 ±0.13 


4.92 ±0.13 


4.98 ±0.13 


5.00 ± 0.06 


Us 

1 


nd 

1 


Full 

3.57 ± 2.22 


optimized binnin 
X^ Bayesian 
3.57 ± 2.22 3.73 ± 2.43 


g 

Poisson 
4.05 ± 2.70 


True PDF 

3.81 ±2.44 


1 


10 


5.05 ± 1.80 


5.29 ±1.72 


3.99 ±1.75 


2.67 ± 1.89 


4.53 ± 0.49 


1 


100 


4.95 ± 1.69 


7.00 ±7.66 


3.47 ±1.65 


2.76 ±1.81 


4.92 ±0.15 


1 


1000 


3.62 ±1.49 


5.86 ± 1.92 


2.30 ±1.31 


2.51 ± 1.61 


4.99 ± 0.07 


10 


1 


3.42 ±2.13 


3.48 ± 2.20 


3.46 ±2.11 


3.31 ±2.18 


3.71 ±2.34 


10 


10 


4.62 ± 0.78 


4.62 ± 0.78 


4.43 ± 0.80 


3.22 ± 1.48 


4.53 ± 0.48 


10 


100 


4.98 ± 0.63 


5.20 ± 0.68 


4.66 ± 0.62 


4.58 ± 0.72 


4.91 ±0.15 


10 


1000 


5.01 ± 0.71 


7.05 ±1.14 


4.26 ± 0.62 


4.28 ± 0.75 


4.99 ± 0.06 


100 


1 


3.48 ±2.16 


3.63 ± 2.35 


3.53 ±2.22 


3.47 ±2.17 


3.90 ±2.52 


100 


10 


4.55 ± 0.56 


4.48 ± 0.63 


4.44 ± 0.57 


4.34 ± 0.98 


4.54 ± 0.50 


100 


100 


4.92 ± 0.27 


4.92 ± 0.27 


4.85 ± 0.27 


4.92 ± 0.26 


4.92 ±0.16 


100 


1000 


4.99 ± 0.24 


5.09 ± 0.25 


4.90 ± 0.24 


5.00 ± 0.25 


4.99 ± 0.06 


1000 


1 


3.33 ±2.12 


3.53 ± 2.39 


3.51 ±2.32 


3.33 ±2.12 


3.82 ± 2.46 


1000 


10 


4.52 ± 0.51 


4.38 ±0.75 


4.42 ± 0.53 


4.49 ± 0.60 


4.53 ± 0.49 


1000 


100 


4.91 ± 0.20 


4.86 ±0.21 


4.87 ±0.19 


4.91 ±0.20 


4.91 ±0.16 


1000 


1000 


4.99 ±0.11 


4.99 ±0.11 


4.97 ±0.11 


5.00 ± 0.11 


4.99 ± 0.06 



Table 1 : Values of /i reconstructed by various methods. Intervals are shown with mean and rms of the reconstructed values 
for 1000 of random data drawings. The result of the first four methods (i.e., except Poisson with true PDF shown in the 
last column) closest to /io = 5 is shown with bold font. For all five methods intervals that do not enclose the true value of 
/io = 5 (i.e., those that are significantly biased) are shown with italic font. The column designated with "Full" corresponds 
to the likelihood description of this paper, section |3] 



10 Conclusion 

In this paper we consider a problem of describing data with limited-statistics simulation sets. After a brief review of the 
usual methods we presented an alternative approach to the likelihood-based comparison of data with multiple simulation 
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sets. In our tests the new approach appears to improve the resuhs of model parameter fits, mainly reducing the bias and 
uncertainty of the fitted parameters. We remark that depending on the statistics of the data and available simulation, and on 
the nature of the experiment other approaches may occasionally produce a better result, so we encourage testing of several 
hkelihood expressions as we did in section|9] 

The method of this paper as described in sections |2]|4] was developed as part of calibration work reported in ^. The 
extension to weighted simulation described in section|5]was proposed for use in several analyses of IceCube, and was shown 
to yield slightly better results (e.g., in [4|: better limits, due to narrower distributions of test statistic). The method of section 
|6]is recently applied to the problem of unfolding energy losses along the muon track in ice with direct re-simulation; the 
correct account of the model error term is important since the model error was not completely removed by our calibration 
efforts (see |3J). 
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